T>aa/ 


{MASA-CB- 177 1 H 3) MODELING CF IfiANSIEtil HEAT 
FIFE GPisRATICk Semiannual Status Report, 19 
Auq. 19£S - l£ let. 1966 (Gectqia Inst, of 
Tech.) 5 2 p CSC! 2QD 


MODELING OF TRANSIENT HEAT PIPE OPERATION 


ii86-3C0S3 

Unclas 

G3/34 43557 


SEMIANNUAL STATUS REPORT 


By 

Gene T. Colwell 
James G. Hartley 


Prepared for 

NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 
LANGLEY RESEARCH CENTER 
HAMPTON. VIRGINIA 23665 


Under 

NASA Grant NAG-1 -392 


June 1 986 


GEORGIA INSTITUTE OF TECHNOLOGY 

A UNIT OF THE UNIVERSITY SYSTEM OF GEORGIA 
SCHOOL OF MECHANICAL ENGINEERING 
ATLANTA, GEORGIA 30332 


NASA GRANT NAG-1 -392 


MODELING OF TRANSIENT HEAT PIPE OPERATION 


By 


Gene T. Colwell and James G. Hartley 
School of Mechanical Engineering 
Atlanta, Georgia 30332 


Submitted to 


National Aeronautics and Space Administration 
Langley Research Center 
Hampton, Virginia 23665 


NASA Technical Officer 
Charles J. Camarda 
Mail Stop 246 


Period Covered 

August 19, 1985 through February 18, 1986 


June 17, 1986 


TABLE OF CONTENTS 


Page 

INTRODUCTION 1 

GOVERNING DIFFERENTIAL EQUATIONS 2 

ANALYTICAL SOLUTIONS... 3 

The Temperature of a Semi -infinite Body 3 

Phase Change of a Semi -infinite Body 4 

DESCRIPTION OF NUMERICAL METHODS 4 

Finite Element Formulation 4 

Time-Stepping Techniques 5 

Latent Heat Evolution Schemes 7 

Boundary Conditions.. 13 

SOME ILLUSTRATIVE TESTS 13 

Example 1. Temperature of a Semi -Ini fi nite Body 13 

Example 2. Solidification of a Semi -Infinite Body in Liquid 15 

Example 3. Phase Change of Sodium 21 

REFERENCES 29 


- ii 



INTRODUCTION 


This report summarizes progress made on NASA Grant NAG-1-392 during the 
period August 19, 1985 to February 18, 1986. The goal of the project is to 
develop mathematical models and associated solution procedures which can be 
used to design heat pipe cooled structures for use on hypersonic vehicles. 
The models should also have the capibility to predict off-design performance 
for a variety of operating conditions. It is expected that the resulting 
models can be used to predict startup behavior of liquid metal heat pipes to 
be used in reentry vehicles, hypersonic aircraft and space nuclear reactors. 

The program work statement includes the following tasks. 

I. To investigate the physics of melting and freezing of the 

working fluid in liquid metal heat pipes, and of other 
performance limitations which are presently not well 

understood, and to develop equations, adequate for design, to 
predict associated performance limits. 

II. To evaluate the Langley finite difference heat pipe design 
computer program and the equations upon which it is based and 
suggest improvements or additions to upgrade the. program. 

III. To work with the aerospace contractor who wins the heat pipe 
study contract to upgrade the contractor's heat pipe analysis 
capability. 

IV. To work with NASA to define the required instrumentation for 
the heat pipe model which will be built under the heat pipe 
study contract and to develop a test plan for the heat pipe. 

Previous status reports in this series (February 1985 and September 1985) 
have covered the development of governing differential equations for the 
startup behavior of initially frozen liquid metal heat pipes. The current 
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report sunmarizes work to date related to numerical solutions of governing 
differential equations for the outer shell and the combination capillary 
structure and working fluid. The working fluid may be solid (frozen), liquid, 
or partly liquid and partly solid. Finite element numerical equations using 
both implicit, explicit and combination methods have been examined. Both 
enthalpy and specific heat approximations to melting have been examined. The 
various numerical solutions have been verified using available analytical 
solutions and experimental data. Work is now in progress aimed at modeling 
vaporization and condensation processes and behavior of the vapor region 
during startup. The next status report will cover this work. 

GOVERNING DIFFERENTIAL EQUATIONS 

The problem under consideration is the initially liquid state at a 
constant temprature T Q , which is greater than the freezing temperature T m . 
For time t > 0, some surfaces are maintained at a constant temperature T,. that 
is lower than T m and the others are insulated. A problem with a phase change 
of a substance from one state to another is mathematically described as 
follows 
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where subscripts s and i indicate the solid and liquid phases, respectively, 
and n represents the unit outward normal direction at the interface. 
Temperatures T., T_, T denote respectively the initial, freezing, and 
boundary tempreaturs. H is the latent heat of fusion, and S is the 
interface postion in the X, Y, t domain. 

m 

ANALYTICAL SOLUTIONS 

The problems of cooling and phase change of semi -infinite body with 
constant boundary temperature has been solved by several authors. In the 
present project, these solutions will serve as reference values for numerical 
calculation. 

The Temperature of a Semi -infinite Body 

The exact solution to this problem is given by Luikou [1] as 


where 


T(x.t) - T 
T o * T w 
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Phase Change of a Semi -infinite Body 


The analytical solution of phase change is given by Ozisik [2] as 


T w erf[x/2(a c t) 1/2 ] 


m w 


erf (x) 
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T q erfc[x/2(a £ t) 1//2 ] 
T o " erfc[x(a s /a Jl ) 1/2 ] 


for solid region 


for liquid region 
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and the location of the interface is 
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where 


, _ S(t) 

' 2(a $ t) 1 / 2 

DESCRIPTION OF NUMERICAL METHODS 
Finite Element Formulation 

The Galerkin weighted residual method is used to derive finite element 
formul ations. Within each element, the unknown function T may be approximated 
at any time t by the relationship 

T(x,y,t) = ^ N i ( x »y) -^(t) (11) 

i=l 

where n is the number of nodes assigned to the element, T- are the descrite 
nodal values of T, and are the shape functions. 

By applying the Galerkin weighted residual method to equations (1) and 
(2), and sujming the contributions from each element and boundaries, assembly 
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of the nonlinear transient element equations can be expressed in matrix form 
as 


[C] { T> + [K] {T> = {F} 


( 12 ) 


where [C] and [K] are, respectively, the capacitance and conductance matrices. 
Time-Stepping Techniques 
Implicit Method 

Let t n denote a typical time in the response so that t n+1 = t n + At 
where At is the time increment, and an intermediate time t m within each time 
step may be expressed as 


t = t + mAt , 0 < m < 1 

m 


Then, Equation (12) at t is given as 


[C] {T} m + [K] {T} m = (F) 


By using a Taylor expansion. 


{T} n = {T} - 4r {T} • (mAt) + High Order Terms 

m o u m 


{T} n+ * = {T} + 4r {Tl • (mAt) + Higher Order Terms 

m dt m 


Subtracting Equation (15) from (16) yields 


<L {T } = - (T) n 

dt * ^m At 


(13) 


(14) 


(15) 

(16) 
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After the higher order terms of Equation (15) are neglected, substitution 


of Equation (15) into (17) gives 

{T) m = (1 - m){T} n + m{T} n+1 (18) 

Substitution of Equation (17) and (18) into the assembly Equation (14) 
yields 


( ^ + «i[K]) (T)" +1 * < - 0 - ")[K])(T}" + {F) (19) 

Since the values of [C], [K] and {F} depend on {T}^ , a choice from among 
the values m = 0 1/2, 2/3 and 1, respectively, yields explicit Euler forward- 
difference, Crank-Nicol son center-difference, Galerkin, and fully implicit 
backward-difference formulations. The fully implicit backward-difference is 
unconditionally stable and predicts a smooth decay and is therefore chosen 
with a Newton-Raphson iteration for the implicit method. 

Explicit Method 

The implicit method is much more stable then the explicit method, but 
requires an iteration within each time step. To avoid- iteration, a three 
level scheme proposed by Lees [3] was used by Comini et al . [4-8] and Morgan 
[9], Oscillations have been observed in certain circumstance, so a slightly 
modified form was used to improve stability. This method has been used 

successfully for phase change problems by several investigators. 

Another three level scheme is referred to as the DuPont three level 
schmem. Hogge [10] demonstrated the overall performance of the DuPont method 
to be superior in accuracy and stability to other time stepping schemes in 
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solving the one-dimensional homogeneous equation. 

Thomas [11] compared several time integration schemes such as the Lees, 
the DuPont scheme, and Crank-Nicol son method. He concluded that the DuPont 
three level scheme was clearly superior to that of Lees in both accuracy and 
stability, and that temperature-dependent terms should be evaluated 
using (T} n+ * instead of {T} n+ ^ . By using the DuPont method. Equation (12) 
can be approximated as 

( I [K] + ){T>n+2 = tiT {T}n+1 + ^ {T >" + {F} (20) 

The equations above allow the explicit evaluation of {T} without 
iteration provided that {T} n+1 and {T} n are known. However, this scheme is 
not self-starting, so {T} n+ * at the first time step may be calculated by using 
the implicit method. 

Latent Heat Evolution Schemes 

The principal difficulties in the analysis of the phase change problem 
are that the variation of the heat capacity is very rapid at the interface as 
shown in Figure 1, the position of the moving interface is not known a priori 
and their shapes may be mul ti -dimensional . Thus, physically realistic 
approximation techniques must be used to overcome the difficulties. 

It is convient to divide them into two groups, based on the choice of 
grid used. In the first group, the moving mesh technique continously tracks 
the location of the interface by deforming the grid system to maintain the 
finest mesh in the vicinity of the critical phase change region. This 
technique may be limited to very simple geometries. 

In the second group, a fixed grid technique can avoid tracking down the 
position of the moving interface, but the interface is generally at an unknown 
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location between nodes. Many different classes of methods are available with 
the fixed grid sysstem. 

The first method uses the enthalpy as a dependent variable along with the 
temperature, and may be referred to as the enthalpy method. Since two 
dependent variables are used, the system of algebaric equations are solved by 
iteration. 

The second method treats the latent heat effect accompanying a change of 
phase in terms of a temperature-dependent specific heat or with the use of an 
enthalpy function. 

These methods avoid the moving interface difficulty so that instead of 
continously tracking down the position of the moving interface, the same 
numerical scheme for conduction heat transafer without phase change is equally 
applicable to the phase change period. Then, the postion of the interface can 
be easily determined by linear interpolation of nodal temperatures. 

Specific Heat Method 

This method approximate the rapid variations of heat capacity over the 
phase change temprature, which is artificially defined, instead of using the 
Dirac Delta function as shown in Figure 1. Thus, a volumetric heat capacity, 
which includes the latent heat of phase change, can be expressed as 


C = 


for T < T - aT 
m 

for T > T + AT 
m 


H c» C c + C 

S£ S l 

2aT 2 


for T - aT < T < T + aT 
m m 


where aT is the phase change temperature interval. 

One advantage of the finite element method is the ability to formulate 
contributions for individual elements before putting them together to 
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Figure 2. Element Model Showing Solid, Liquid, 
and Fusion Phase. 
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represent the entire domain. The individual element characteristics are 
determined by the shape functions and nodal values. Thus, the heat capacity 
of the fusion element can be evaluated correctly instead of the direct use of 
the averaged element temperature. 

Consider the phase change element in Figure 2 and assume that the fusion 
phase is within the element. If the variation of the density of the working 
substance around the freezing temperature is quite small, the mass fractions 
can be replaced by area fractions M , M s> for liquid, solid, and fusion 
phases with unit thickness. Otherwise, mass fractions which are calculated 
based on the area and density of each phase is employed. 

For a linear element, the phases within the element are seperated by the 
isothermal lines, T m + aT and T^ - aT , passing through the element. The 
intersecting points of the isothermal lines at the boundaries of the element 
and nodal points are the vertices of the area occupied by each phase. Since 
three-node, linear triangular elements are used, at least one of the phases 
occupies an area having a triangular shape. Then, the area of the traingular 
is obtained by using the known coordinate values of the vertices of the 
triangle. For example, the area fracton of the element as shown in Figure 2 
are expressed as 



M s = 1 - <\ - V 
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where b is the total area of the element with the sign chosen so that the area 
is nonnegative. This new technique, which uses the mass fractions 
Mj. for liquid, solid, and fusion phases, is introduced to correct 
misinterpretation of the specific heat without a vry fine element grid and/or 
time step. 

Enthalpy Function 

When temperature approaches the phase ■ change temperature, the heat 
capacity tends to the Dirac delta function and can not be satisfactorily 
represented across the peak by any smooth function. In Comini et al [4], a 
technique based on the integral of heat capacity with respect to temperature 
is proposed 

H = f T P C dT 

oo 

This is a smooth function of temperature in the phase change zone. 
Therefore, the enthalpy is interpolated rather than the heat capacity in a 
element as following 

n 

H = L N. (X,Y)H. (t) (22) 

i = l 1 1 

where again are the shape functions and are the enthalpy values at nodal 
points. 

From definition, the heat capacity can be expressed as 

pc = dH/dT (23) 
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Thus, the values of the heat capacity can be approximated by evaluating 
the gradient of enthalpy with respect to temperature. Defining the direction 
S to be normal to the interface line, Equation (23) is expressed as 


_ / .3H . _3T . 
pc ‘ ^ as 1 as ) 


where 


and 


= ( M * + — SL )/{ — ) 

K 3x *sx 3Y Vy' M 3S ' 


o = U / II o _ il / il 
sx 3x ' 3s * sy 3y ' 3s * 


il - r f il i 2 + ^ il \h l,z 

3S “ L ' 3X ' V ; J 


(24) 


ay 


Hence, for the entire element, the final expression of the heat capacity 
is given as 


r 1 _ r d . H 3T 3H 3T -i/rr 9T .2 , 3T .2-1 

^ p ^ ~ ^ 3x * 3x 3y * 3y 3x ^ ^ 3y ^ ^ 


(25) 


where, for linear temperature, tringular element. 


U 1 

3x = 2 a (b l H 1 + b 2 H 2 + b 3 H 3^ 


gy = 2 A ( C l H l + C 2 H 2 + C 3 H 3^ 


3X = 2a ( b l T l + b 2 T 2 + b 3 T 3^ 


3Y = 2A ^ C 1 T 1 + C 2 T 2 + C 3 T 3^ 
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The appropriate pc value to be used in the element can be evaluated by 
Equation (25). This expression is developed by Del-Giudice [12]. 

Boundary Conditions 

Since specified temperature boundary conditions are abruptly imposed on 
the boundary surface, the temperature gradient at this region is infinite 
which may cause unstable conditions. This difficulty may be elimilited by 
introducing a layer of fictitious elements with negligible thermal storage 
capacity and very high thermal conductivity, such that the temperature 
gradient at this region artificially becomes of finite value. For this 
purpose, numerical values of the thermophysical properties of the fictitious 
layer are chosen as following: 

pc = 1.26 x 10 4 J/m 3 K 

k = 4.19 x 10 3 W/m K 

SOME ILLUSTRATIVE TESTS 

Example 1. Temperature of a Semi -Ini finite Body 

This example involves pure conduction of heat, without phase change, over 
a semi-infi nite body. It is solved as a two dimensional problem shown in 
Figure 3. This problem is used to test the general capability of the program 
written in FORTRAN V. 

Adiabatic boundary conditions are assumed throughout but at the face x = 
0 the specified boundary temperature (253.0 K) is maintained constant during 
the entire heat transfer process while the initial solid temperature (283.0 K) 
is uniform. The conductivity and specific heat are assigned values of 
al umi num. 
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Figure 3. Finite-Element mesh for transient conduction heat transfer 
problem. 


Both the fully implicit backward-difference and the explicit time 
stepping schemes are tested. The explicit scheme causes small oscillations at 
early time steps. Hence an implicit scheme is employed to eliminate the 
oscillations in the first five time steps. A time step of 10 seconds is used, 
and numerical calculations are terminated as asymtotic conditions are 
approached. 

Numerical results are compared in Figure 4 to the well known analytical 
solution. Both results from the implicit and explicit schemes agree well with 
the analytical solutions. 

Example 2. Solidification of a Semi -Infinite Body in Liquid 

This example as represented in Figure 5 tests the ability of the 
numerical methods to handle the latent heat of phase change which accompanies 
the discontinuity of the properties at the phase change region. 

The same initial temperature and boundary conditions as example 1 are 
used. The initial temperature is higher than freezing (273.0 K). The 
properties of water are assigned to the specific heat and conductivities of 
the solid and the liquid state. The phase change interval, 2 aT is assumed to 
be .5 °K and a time step of 400 seconds is used for Figure 6, 7 and 8. 

Numerical results using the specific heat method and the enthalpy 
function together with a implicit or explicit time stepping scheme are 
compared with the results from an analytical solution [2] in Figures 6 and 7, 
respectively. As shown in Figures 6 and 7, the temperature predictions of the 
solid region by both methods agree well with analytical solutions. In the 
liquid region, the maximum deviation is about 2°K. 

Numerical results predicted by the explicit method using an implicit 
scheme for the initial five time steps are closer to the analytical ones than 
that by the implicit scheme. Also, the explicit method consumes less 
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computational time since iteration is not involved. 

In Figure 8, numerical results calculated by using specific heat and 

enthalpy function with an explicit time- stepping scheme are compred with 
analytical ones. At early time steps, the specific heat method predicts 
better temprature distributions than that of the enthalpy function but for the 
later time steps, numerical results predicted by using the enthalpy function 
gives better agreement to analytical solutions. This method requires much 
less computational time than the specific heat method. In the solid region, 
both methods yield quite good results. 

To predict temperature distributions for longer time, the enthalpy 

function with an explicit time scheme are chosen because the results shown in 
Figures 6, 7 and 8. The same initial and boundary conditions are used, but 
the domain is extended as shown in Figure 5 and a time step of 500 seconds is 
used. 

The calculated results are compared with analytical solutions in Figure 
9. Good prediction is achieved except for the phase change region, but even 
in this region the maximum deviation is less than 2 K. In Figure 10, the 

position of interfaces calculated by linear interpolation of adjacent nodal 

temperatures of the freezing temperature gives good agreement with analytical 
sol utions. 

Example 3. Phase Change of Sodium 

In this example, sodium is used as the phase change material for one and 
two-dimensional regions. The initial temperature of the sodium is 393. K which 
is higher than the freezing temperature (371. K). First kind boundary 
conditions (293. K) are imposed on boundary surfaces. The phase change 
interval, 2&T is assumed to be 2.K and a time step of 5 seconds is employed. 
The finite element meshes are shown in Figures 5 and 11. 
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Figure 10. Position of solidification front at different time 
values referred to in Figure 9. 



Adiabatic Boundary 
Condition 



Figure 11. Two-dimensional mesh used to represent a 
corner region; 242 elements, 144 nodes. 




The enthalpy functions with an explicit time stepping scheme is used for 
the phase change of sodium. In Figures 12 and 13 the temperature 

distributions for pure conduction and phase change in a one-dimension region 
show excellent agreement with analytical solutions. 

For the two-dimensional case, the position of the interface is compared 
with an analytical solution given by Jiji [13]. This analytical solution 

assumes that the ratio of diffusivity for the solid to liquid is unity and 
that the interface at points beyond three times the one-dimensional stationary 
interface in x and y directions is the same as the one-dimensional stationary 
interface position. 

Since the ratio of di ffusi vities is not unity, the numerical solutions 
can not be directly compared with the analytical solutions. However, as shown 
in Figure 14, the solutions are close to each other for early time steps 

because adiabatic boundary conditions at surfaces which do not have first kind 
boundary conditions are effective for both solutions. The difference of 
diffusivities relatively small. As time increases the positions of the 
interface given by numerical calculation advance further than the analytical 
ones and interface lines obtained from the analytical method are not 

perpendicular to boundary surfaces. This means that numerical calculations 
matches adiabatic boundary conditions, but the analytical solution does not. 
Also, the effect of the difference of diffusivities is not negligible since 
quite a lot of mass has changed from liquid to solid. Even though numerical 
solutions do not exactly match analytical solutions because of the different 
conditions of the problem, general trends of the numerical results are 

veri fied . 
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POSITION OF INTERFACE AT DIFFERENT TIMES 


Figure 14. Solidification of sodium in a corner region 
position of interface at different time's, 
a time step of 5 seconds is used. 
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